
cd('.')
addpath('matlab functions')

D0 = readtable('NC Data/SCC_4to5.csv');

OUT = struct();

for sbj = {'math','reading'}

    vamdl = load(sprintf('output/main/vam_%s_matching.mat',sbj{1}));

    tchcoef = vamdl.tchcoef;

    tchcoef_split = cellfun(@(x) strsplit(x,'X'),tchcoef,'unif',0);
    tchcoef_split{strcmp(tchcoef,sprintf('lag_%sscore_sq',sbj{1}))} = repmat({sprintf('lag_%sscore',sbj{1})},[1 2]);
 
    D = D0(:,unique(cat(2,tchcoef_split{:})));

    numpts = 501;
    OUT.(sbj{1}).score = struct('x',linspace(-2.5,2.5,numpts)','y',nan([numpts 1]));
    for i = 1:numpts
        Z = D;
        Z{:,sprintf('lag_%sscore',sbj{1})} = OUT.(sbj{1}).score.x(i);
        Z = cellfun(@(x) prod(Z{:,x},2),tchcoef_split,'unif',0);
        Z = cat(2,Z{:});
        OUT.(sbj{1}).score.y(i) = mean(sqrt(sum((Z*vamdl.VAmomentEstimates.cov).*Z,2)));
    end

    for f = {'female','race_black','poverty','lep'}
        OUT.(sbj{1}).(f{1}) = nan([1 2]);
        for i = 0:1
            Z = D;
            Z{:,f{1}} = i;
            Z = cellfun(@(x) prod(Z{:,x},2),tchcoef_split,'unif',0);
            Z = cat(2,Z{:});
            OUT.(sbj{1}).(f{1})(i+1) = mean(sqrt(sum((Z*vamdl.VAmomentEstimates.cov).*Z,2)));
        end
    end

    g_vals = {[0 0] [0 1] [1 0] [1 1]};
    for g = 1:numel(g_vals)
        Z = D;
        Z{:,'female'} = g_vals{g}(1);
        Z{:,'race_black'} = g_vals{g}(2);
        Z = cellfun(@(x) prod(Z{:,x},2),tchcoef_split,'unif',0);
        Z = cat(2,Z{:});
        OUT.(sbj{1}).race_gender.cond(g) = mean(sqrt(sum((Z*vamdl.VAmomentEstimates.cov).*Z,2)));

        Z = D(D.female==g_vals{g}(1) & D.race_black==g_vals{g}(2),:);
        Z = cellfun(@(x) prod(Z{:,x},2),tchcoef_split,'unif',0);
        Z = cat(2,Z{:});
        OUT.(sbj{1}).race_gender.uncond(g) = mean(sqrt(sum((Z*vamdl.VAmomentEstimates.cov).*Z,2)));
    end

end

plotformat = {'interpreter', 'latex','FontWeight','normal','FontSize',10,'FontName','Times'};
axisformat = {'FontUnits','points','FontWeight','normal','FontSize',10,'FontName','Times'};

%-- Plots variance by lag score --
FOOTER = .08;
OUTPOS = subplot_outpos([1 1],0,FOOTER,.02);
FIGRATIO = 3/4;

close all
pos = get(gcf,'Position');
pos(3) = sqrt((pos(3)*pos(4)/FIGRATIO));
pos(4) = pos(3)*FIGRATIO;
set(gcf,'Position',pos)
    
hold on
plot(OUT.math.score.x,OUT.math.score.y,'-k')
plot(OUT.reading.score.x,OUT.reading.score.y,'--k')
xlabel('Lag Score',plotformat{:})  
set(gca,'YLim',[.1 .35],'YTick',(.1:.05:.35),'XLim',[-2.5 2.5],'XTick',-2.5:.5:2.5,'OuterPosition',OUTPOS{1,1},axisformat{:})
hL = legend({'Math','Reading'},'Orientation','horizontal','location','southoutside',plotformat{:});
set(hL,'Position',[(1-hL.Position(3))/2 (FOOTER-hL.Position(4))/2 hL.Position(3) hL.Position(4)]);

set(gcf,'PaperPositionMode', 'manual','PaperUnits','inches'...
        ,'PaperPosition',[0 0 [1 FIGRATIO]*4.5],'PaperSize',[1 FIGRATIO]*4.5)
print('-dpdf','tables and figures/graphics/va_sd_lagscore.pdf','-r2000')
close all  

%-- Bar Chart Demographics --
FOOTER = .1;
DIM = [1 2];
OUTPOS = subplot_outpos(DIM,0,FOOTER,.02);
FIGRATIO = .55;

close all
pos = get(gcf,'Position');
pos(3) = sqrt((pos(3)*pos(4)/FIGRATIO));
pos(4) = pos(3)*FIGRATIO;
set(gcf,'Position',pos)

i = 1;
for sbj = {'math','reading'}
    X = [OUT.(sbj{1}).female; OUT.(sbj{1}).race_black; OUT.(sbj{1}).poverty; OUT.(sbj{1}).lep];
    X = flip(X,2);

    [COL,ROW] = ind2sub(flip(DIM),i);
    subplot(DIM(1),DIM(2),i)
    b = bar(X);
    b(1).FaceColor = [.1 .1 .1];
    b(2).FaceColor = [.9 .9 .9];
    title([upper(sbj{1}(1)) sbj{1}(2:end)],plotformat{:});
    xticklabels({'Female','Black','FRL','LEP'})
    if i==1
        YLIM = [0.17 .27];
    elseif i==2
        YLIM = [0.1 .2];
    end    
    set(gca,'box','off','YLim',YLIM,'OuterPosition',OUTPOS{ROW,COL},axisformat{:});
    i = i + 1;
end
hL = legend({'Yes (1)','No (0)'},'Orientation','horizontal',plotformat{:});
set(hL,'Position',[(1-hL.Position(3))/2 (FOOTER-hL.Position(4))/2 hL.Position(3) hL.Position(4)]);

set(gcf,'PaperPositionMode', 'manual','PaperUnits','inches'...
        ,'PaperPosition',[0 0 [1 FIGRATIO]*4.5],'PaperSize',[1 FIGRATIO]*4.5)
print('-dpdf','tables and figures/graphics/va_sd_demo.pdf','-r2000')
close all  


%-- Bar Chart Race Gender --
FOOTER = .16;
DIM = [1 2];
OUTPOS = subplot_outpos(DIM,0,FOOTER,.02);
FIGRATIO = .65;

close all
pos = get(gcf,'Position');
pos(3) = sqrt((pos(3)*pos(4)/FIGRATIO));
pos(4) = pos(3)*FIGRATIO;
set(gcf,'Position',pos)

i = 1;
for sbj = {'math','reading'}
    X = [OUT.(sbj{1}).race_gender.cond; OUT.(sbj{1}).race_gender.uncond];

    [COL,ROW] = ind2sub(flip(DIM),i);
    subplot(DIM(1),DIM(2),i)
    b = bar(X);
    b(1).FaceColor = [.1 .1 .1];
    b(2).FaceColor = [11/30 11/30 11/30];
    b(3).FaceColor = [19/30 19/30 19/30];
    b(4).FaceColor = [.95 .95 .95];
    title([upper(sbj{1}(1)) sbj{1}(2:end)],plotformat{:});
    xticklabels({'',''})
    if i==1
        YLIM = [0.17 .27];
    elseif i==2
        YLIM = [0.1 .2];
    end    
    set(gca,'box','off','YLim',YLIM,'OuterPosition',OUTPOS{ROW,COL},axisformat{:});
    x = text(1:2,(min(ylim)-0.0095)*ones(1,2),{{'Other';'Variables Fixed'} {'Unconditional'}},'HorizontalAlignment','center',axisformat{:});
    i = i + 1;
end
hL = legend({'Non-Black Male','Black Male','Non-Black Female','Black Female'},'NumColumns',2,'Orientation','horizontal',plotformat{:});
set(hL,'Position',[(1-hL.Position(3))/2 (FOOTER-hL.Position(4))/2 hL.Position(3) hL.Position(4)]);

set(gcf,'PaperPositionMode', 'manual','PaperUnits','inches'...
        ,'PaperPosition',[0 0 [1 FIGRATIO]*4.5],'PaperSize',[1 FIGRATIO]*4.5)
print('-dpdf','tables and figures/graphics/va_sd_racegender.pdf','-r2000')
close all  



